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IMPACT AND PENETRATION SIMULATIONS 
FOR COMPOSITE WING-LIKE STRUCTURES 


— Final Report for NAG-1-1858 - 

by 

Norman F. Knight, Jr. 

Department of Aerospace Engineering 
College of Engineering and Technology 
Old Dominion University 

SUMMARY 

The goal of this research project was to develop methodologies for the analysis of 
wing-like structures subjected to impact loadings. Low-speed impact causing either no 
damage or only minimal damage and high-speed impact causing severe laminate damage 
and possible penetration of the structure were to be considered during this research effort. 
To address this goal, an assessment of current analytical tools for impact analysis was 
performed. Assessment of the analytical tools for impact and penetration simulations with 
regard to accuracy, modeling, and damage modeling was considered as well as robustness, 
efficient, and usage in a wing design environment. Following a qualitative assessment, 
selected quantitative evaluations will be performed using the leading simulation tools. 
Based on this assessment, future research thrusts for impact and penetration simulation 
of composite wing-like structures were identified. 

BACKGROUND 

The design of aerospace structures generally results in lightweight structural de- 
signs which exploit advanced structural materials and fabrication concepts. The goals 
of aerospace structural design are to meet design requirements based on the operating 
conditions and flight envelope of the vehicle, service life and damage tolerance, and man- 
ufacturing cost. Over their intended life cycles, aerospace structures may be subject to 
internal pressure loads, thermal cycling, bending, axial, and shear loads, impact, and fa- 
tigue. Aerospace structures frequently involve flat and curved, stiffened and unstiffened 
panels, with and without cutouts, that are interconnected by frames, stringers and bulk- 
heads. In addition to satisfying the requirements of the normal operating environment, 
the design should also be analyzed to assess the ability of the structure to contain dam- 
age due to impact and penetration, thus establishing, in part, the residual strength and 
crashworthiness of the vehicle. Improved crashworthiness will increase the probability of 
survival for the passengers and crew. 

Many researchers have experimentally examined the susceptibility of selected com- 
posite laminates to low-speed impact damage using relatively small-scale specimens. Oth- 
ers are investigating the dynamic crushing response of composite structures for improved 
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crashworthiness. In addition, researchers have been taking an analytical approach to study- 
ing the impact response of composite structures using a Rayleigh-Ritz approach, a finite 
element approach, or an integral equations approach. Representative papers in these areas 
are available in Reference 1. 

The solution approach commonly used for this class of problems is the solution to 
the nonlinear transient dynamic analysis problem. This approach requires the solution of 
the equations of motion as a function of time during the impact/penetration and/or crash 
event. Various spatial approximations or discretizations axe used giving rise to a semi- 
discrete system of equations (i.e., ordinary differential equations in the time domain). The 
conventional solution methods for direct time integration are the implicit time integration 
methods and the explicit time integration methods. Each solution method is described in 
the subsequent sections. Evaluation of high-performance equation solvers and a comparison 
of implicit and explicit procedures are given in References 2 and 3, respectively. 

Conventional Implicit Solution Methods 

The semi-discrete finite element equations of motion for time increment n + 1 may be 
written as 


MD n+1 + CD n+1 + F(D n+1 ) = P n+1 (1) 

where M is the mass matrix, C is the damping matrix, F is the internal force vector, and P 
is a vector of external loads. The parameters D, D, and D represent acceleration, velocity 
and displacement vectors, respectively. It should be noted that internal force vector F is 
a function of the displacements as 

F(D n+1 ) = K 0 D n+1 + Q(D n+1 ) (2) 

where K 0 represents the linear stiffness matrix and Q is a vector of nonlinear terms 
which may arise from geometric/material nonlinearities, contact /boundary conditions, and 
follower forces. The solution of equation (1) may be obtained by employing any one of 
many direct time integration algorithms. A popular choice is the implicit Newmark method 
which provides the following approximations for the velocity D n+1 and acceleration D n+1 
vectors at time step n + 1: 

6 ” +1 = Zs< D " +1 - D ”> - (f - 0 6 " - ($ - 0 hi> " (*o 

6 " +i = ^ 

where h represents the time step (i.e., h = t n +\ — f n ), and the parameters o and (3 govern 
properties of the algorithm. Unconditional stability is achieved when a = | and P = 7 
which results in the constant average- acceleration version of the Newmark method. Using 
these definitions for ot and P and substituting equations (3) into equation (1), the equations 
of motion may now be expressed as 
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The accompanying equations for F m n and F c n which represent historical inertial and 
viscous force vectors from the previous time step, respectively, are 


Fra" = M 


4 4 • 

^D»+ -D"+D" 


(5a) 


F c n = C Qd" + D" j (5b) 

The right-hand side of equations (4) is now independent of displacements at time n + l. For 
most complex applications such as an impact/penetration or crash analysis, these equations 
are nonlinear owing to equation (2) for the internal forces. As such, the Newton- Raphson 
iterative technique is generally applied at each time step, and the desired displacement 
solution D n+1 is obtained through incremental updates for the kth iteration as 


d;+; = dj +1 + ad; +1 ( 6) 

where and denote the current (kth) and updated ( k + 1st) iterative displace- 

ment estimates at time n + l for the Newton- Raphson iteration, respectively. The dis- 
placement increment AD£ +1 is obtained by solving a system of linearized equations which 

must be evaluated, assembled, factored or decomposed, and solved for each iteration. That 
is, 


K£ +1 AD£ +1 = R£ +I (7) 

where the “effective” tangent stiffness matrix K and the residual or out-of-balance force 
vector R are given by 


KjT 1 



+ h M + l c 


( 8 ) 


R n+1 = pn+1 + Fm n + _ p^n+l) (9) 


The computational intensity of nonlinear transient dynamic analyses using an implicit 
time integration method derives from the fact that each iteration requires the assembly 
and factorization of K, and that this procedure must be repeated for each time step. Use 
of a modified Newton-Raphson scheme may reduce the number of factorizations; however, 
for highly nonlinear applications such as an impact /penetration or crash analysis, the asso- 
ciated gains are generally outweighed by a substantial increase in the number of iterations. 
It should be noted that the development just outlined is based on the more general case of 
transient dynamic analysis. For static analysis, the time dependence would be removed, 
and the inertia and damping terms in equation (1) are neglected. Solutions for the static 
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analysis also require the Newton-Raphson iterative technique and are also computationally 
intensive. 

The implicit scheme considered herein is unconditionally stable, and the time step 
size is limited only by accuracy considerations. This generally permits the use of a large 
time step. However, several Newton-Raphson increments are often required for each time 
step, and this number can increase substantially for problems which are highly nonlinear or 
exhibit structural instabilities. Each Newton-Raphson iteration requires the solution of a 
linearized system of equations of the form Ax = b, where A, the effective tangent stiffness 
matrix, is sparse, symmetric, and may or may not be positive definite. The solution is 
obtained using either direct or iterative methods for solving systems of linear algebraic 
equations. 

Direct solvers are based on elimination or factorization techniques and are variants 
of Gaussian elimination. The most widely used and heavily researched technique is the 
Cholesky method. These techniques are computationally intensive (i.e., 0(n 3 ) floating- 
point operations for the factorization step alone) but reliable - a solution is guaranteed 
after execution of a fixed number of operations. However, they are difficult to parallelize 
efficiently due to the sparsity of the system of equations, and the communication needs. 
Furthermore, they are very memory intensive since the full system of equations is required. 

Iterative solvers for systems of linear algebraic equations axe successive approximation 
techniques based on an initial guess for the solution. One such method which is being heav- 
ily researched is the Preconditioned Conjugate Gradient or PCG (e.g., see Reference 2). 
This method may be implemented on an element-by-element basis resulting in low mem- 
ory and communication requirements, and is therefore well-suited for parallel processing. 
However, its reliability is sensitive to the conditioning of A. Ill-conditioning due to geo- 
metric and material disparities can lead to extremely slow convergence, no convergence, 
or even convergence to the wrong results. 

In summary, the kernel of the implicit approach is the repeated evaluation, assembly, 
and solution of a linearized system of equations. Conventional direct solvers are fast and 
reliable on sequential computers; however, their communication and memory requirements 
lead to poor scalability in the context of parallel processing. Conventional iterative solvers 
are very amenable to parallel processing; however, they suffer from convergence problems 
due to ill-conditioning. Thus, the need exists for a reliable alternative approach which is 
inherently suited for parallel processing. This alternative approach should accommodate 
the computational requirements associated with very difficult and complex structures such 
as composite flight-vehicle structures. It should also exploit the capabilities of a wide range 
of high-performance parallel computer architectures. 

Conventional Explicit Solution Methods 

The use of an explicit, rather than implicit, time integration technique to solve the 
semi-discrete system of equations given by equation (1), is attractive for parallel computa- 
tions. Central-difference approximations are typically employed for the temporal deriva- 
tives. That is, 
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(10a) 


D" +1 = J-(D n+1 - D"” 1 ) 

2 h 

D n+1 = ^-(D n+1 - 2D" + D n_1 ) (106) 

Substituting these approximations into equation (1), the following equation for displace- 
ments at time n + 1 may be obtained 

0^ M+ ^ C ) Dn+1 = pn_ F ( Dn ) + /| MDn - ^ MD n_1 + ^CD"- 1 (11) 

where h is the time step. These equations are linear (even for nonlinear problems), and the 
left-hand-side matrix is constant unless the time step h changes during the solution process. 
Also, if diagonal mass and damping matrices are used, these equations represent an un- 
coupled system of algebraic equations in which each solution component may be computed 
independently. For transient dynamic analysis, a time history of displacements (system 
response) is sought. Mass and damping vectors which best model the physical properties 
of the system are used. Techniques for estimating the maximum allowable time step axe 
available, such that the time step size may change during the transient dynamic analysis. 
As such, explicit time integration techniques axe attractive candidates for implementation 
on parallel computers. These techniques generally have low memory and communication 
requirements but are also only conditionally stable numerically. Effective solution to the 
static problem on parallel computer systems still remains the same; however, recent work 
on adaptive dynamic relaxation procedures is very promising (e.g., see Reference 4). 

OVERALL RESEARCH GOALS 

The goal of this research activity was to develop and assess methodologies for the sim- 
ulation of impact and penetration of composite wing-like structures. This effort focussed 
on low-speed, single-site impact damage with the overall goal of being able to simulate 
general impact damage at multiple sites for large built-up composite structures. Four 
primary objectives were originally included in this overall research activity: 

1. Identification of fundamental mechanics issues associated with impact and pene- 
tration simulations for composite structures. 

2. Identification of analytical tools for these simulations and assessment of their 
capabilities. 

3. Develop an intelligent computational system for modeling and analysis of wing- 
like structures subjected to impact and penetration. 

4. Application of these methodologies to the analysis of composite wing-like struc- 
tures. 

Since this research activity was originally proposed for multiple years of work and 
due to program re-direction only lasted for one year, only the first two objectives were 
addressed. 
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RESEARCH OBJECTIVES AND FINDINGS 


It was anticipated that this research project would involve several years of research 
that would result in an accurate simulation tool for composite structure design for impact 
loading. The goal of this work was to establish the state-of-the-art in simulation tools for 
impact analysis of composite structures. Benchmark cases were used to establish compara- 
tive assessments and limited comparisons with existing test data were performed. Specific 
research objectives for this grant were as followed: 

• Identify key modeling and analysis issues needed to simulate the physics of impact 
and penetration of composite wing-like structures. 

• Assess selected leading large deformation impact and penetration simulation tools for 
applicability to composite wing structures. 

• Perform parametric studies using the leading simulation tools to identify key modeling 
aspects and appropriate analysis options. 

• Identify key research thrusts for developing advanced techniques for impact and pen- 
etration simulations. 

The principle investigator and two graduate students were involved in this research. 
One graduate student completed his Master’s thesis (see Reference 5) using DYNA3D 
as part of this effort, and a second student initiated some follow-on activities using LS- 
DYNA3D that contributed to identifying research thrusts for future work. A copy of the 
Master’s thesis has been sent to the grant monitor under a separate cover letter. This 
report summarizes those findings and those research thrusts identified in the follow-on 
effort. 

Identify Key Modeling and Analysis Issues 

To simulate the impact or penetration event accurately, the simulation tool must in- 
corporate several key features. First, the analysis should account for large deformation 
effects and possibly tearing of the structure. Second, the simulation tool should include 
capabilities for modeling both 2-D and 3-D geometries, laminated structures with and 
without damage, contact, frictional interfaces, and mesh refinement strategies. The com- 
posite damage modeling is dependent on the failure modes and mechanisms implemented 
in the analysis tool. New failure modes and damage models may need to be incorporated 
as they are developed. The contact algorithm should incorporate sliding frictional inter- 
faces as well as possible surface separation or complete penetration. Third, the analysis 
most likely will involve an explicit direct-time integration procedure; however, an optional 
implicit direct-time integration procedure should be available for long duration response 
predictions beyond the initial impact event. Only a preliminary exploratory study of the 
penetration mechanics problem was done. The following issues were identified: 

• Large deformation effects must be included in the simulation. The results for com- 
posites may still involve small strains but the kinematics involved are large deflections 
and rotations. 

• Impact or crash events of metal structures typically result in “folding” of thin sheets 
or large strain effects leading to material failures. In composites, such behavior results 
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in delaminations, cracks and brittle failures. Tearing of material requires a modeling 
and analysis capability much different than a standard nonlinear transient analysis 
simulation. 

• Erosion of elements or tearing of connected elements must be a feature of the analysis 
tool in order to simulate the progression of through-the-thickness damage. Two of the 
modeling approaches are indicated in Figure 1. Some methods conserve mass (e.g., 
node release) but require a priori knowledge and modeling of the damage growth 
directions. Other methods (e.g., element erosion and elimination) use material failure 
modes to determine damage growth directions but require localized refinement in order 
to minimize the mass loss due to element erosion. 

• Modeling features need to include both a 2-D and 3-D capability depending on the 
fidelity of the response prediction desired. In the 2-D plate/shell simulations, the 
through-the-thickness damage is difficult and computationally intensive to simulate. 
Since these formulations generally assume inextensibility in the thickness direction, 
interlaminar stresses are not predicted even though they may be drivers in the failure 
mode. In the 3-D solid simulations, the through-the-thickness damage can be modeled 
but the spatial discretization in the thickness direction has a direct impact on the 
planar spatial discretization (i.e., element aspect ratio). The interlaminar normal 
stresses play a dominate role in damage initiation and growth for laminated composite 
structures and sandwich structures with soft core materials. 

• Existing failure models in most of the simulation tools include extensive elasto-plastic 
material models for metal structures and only point-stress failure analyses with ply 
discounting for laminated composite structures. A capability to add user-defined 
material models is needed in order to implement and assess evolving failure models for 
laminated composites and sandwich structures. One such study of progressive failure 
analysis methods for laminate composite structures under static loading conditions is 
given in Reference 6. The need for reliable material data, a failure model and material 
degradation models is clearly established. 

• Simulation tools should provide both explicit and implicit direct time integration 
methods. The initial impact event and the response immediately after impact, includ- 
ing possible penetration, occur over a very short time interval and is best analyzed 
using an explicit solution method. However, for the long time response prediction, an 
implicit solution method is desirable in order to move forward in time. 

Assess Impact and Penetration Simulation Tools 

A review of several simulation tools was performed. Several of the leading commer- 
cially available finite element codes were included in this assessment including DYNA3D, 
LS-DYNA3D, MSC/DYTRAN, and STAGS. ODU has the LLNL version of DYNA3D and 
NIKE3D and is a member of the LLNL Collaborator Program. LS-DYNA3D is available 
on one Unix workstation at NASA Langley, and ODU obtained an academic license of 
LS-DYNA3D to support this work. ODU negotiated with MSC to extend our existing 
license for NASTRAN and PATRAN to include MSC/DYTRAN, but it was not possible 
to obtain and install the software during the grant period. MSC/DYTRAN is available on 
a limited basis from a single computer at NASA Langley. STAGS is available on the NASA 
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computers, and ODU did not obtain a copy from COSMIC for our on-site computational 
systems primarily due to the limited time period for the grant. Most of these codes, if 
not all, already have several of the analysis capabilities needed for impact/penetration or 
crash analysis of metal structures. For composites, the damage models and treatment of 
historical data associated with nonlinear material response is the most critical. Current 
damage models are based on the observed failure modes which axe not documented for 
large-scale composite built-up structures. Developing and/or incorporating new material 
models will most likely be required in any of the simulation tools. 

Currently, several analysis codes are available from commercial companies or gov- 
ernment laboratories for simulating impact, crash, and penetration events. For the most 
part, these analysis codes focus on the impact event for metal structures in order to assess 
and improve the crashworthiness of automobiles and aircraft, and on ballistic impacts for 
military applications. Over the past 15 years, several comparisons have been made to 
document their features, capabilities, and performance (e.g., see References 7-11). The 
leading analysis tools for such simulations are DYNA3D, LS-DYNA3D, NIKE3D, ANSYS, 
ADINA, WHAMS, DYCAST, MSC/DYTRAN, ABAQUS/Explicit, and STAGS. A brief 
description of each code is given herein. 

• DYNA3D - This code was developed by Lawrence Livermore National Laboratory 
(LLNL) and represents a fully-explicit, nonlinear, transient dynamic, finite-element 
code that models large deformations of nonlinear materials and contact. Sliding con- 
tact with friction and voids (i.e. , separation after contact) is permitted. The integra- 
tion is based on an explicit central- difference time integrator with automatic time step 
adjustment for efficiency and numerical stability. Very limited support for this code 
is available from LLNL. 

• LS-DYNA3D - This code is available from Livermore Software Technology Corporation 
and represents the commercial version of DYNA3D from LLNL. Enhanced analysis 
features and user interfaces are provided. Impact, penetration, crash, and airbag 
deployment are all possible analysis options in this code. 

• NIKE3D - This code was developed by Lawrence Livermore National Laboratory 
and represents a fully-implicit, nonlinear, transient dynamic, finite element code that 
models large deformations of nonlinear materials and contact. Sliding contact with 
friction and voids (i.e., separation after contact) is permitted. The integration is based 
on an implicit time integrator with automatic time step adjustment for efficiency and 
accuracy. 

• MSC/DYTRAN - This code is available from the MacNeal Schwendler Corporation as 
a general-purpose, nonlinear, transient dynamic, finite-element code which includes 
both a Lagrangian and an Eulerian formulation. This code shares its origins with 
the DYNA3D code. MSC/DYTRAN also offers a fluid and fluid/structure interaction 
analysis feature based on PISCES. Impact, penetration, crash, and airbag deployment 
axe all possible analysis options in this code. 

• ABAQUS/Explicit - This code is available from Hibbitt, Kaxlsson, and Sorenson, Inc. 
(HKS) and represents a fully-explicit, nonlinear, transient dynamic, finite-element 
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code that models nonlinear materials and contact. 

• ANSYS - This code is available from ANSYS, Inc. (formerly Swanson Analysis Sys- 
tems, Inc.). It is a general-purpose, finite-element code for linear and nonlinear, 2-D 
and 3-D structures. 

• ADINA - This code is available from ADINA R&D, Inc. It is a general-purpose, 
finite-element code for linear and nonlinear 2-D and 3-D structures. It claims that 
the overall reliability, efficiency, and accuracy of ADINA for state-of-the-art practical 
analysis distinguishes it from other finite-element codes. The code and company were 
developed by Professor K. J. Bathe of MIT and his associates. 

• DYCAST - This code was developed by Northrop-Grumman Corporation under NASA 
sponsorship. It contains both explicit and implicit direct time integration methods 
and is focused on crashworthiness predictions for aircraft, helicopter and automotive 
applications. This code has been developed over a number of years with Dr. Allan 
Pifko as the main technical developer. The capabilities axe heavily focused on crash 
simulation. However, due to limited development manpower, the code lacks many 
features needed for general crash and impact simulations of composite structures, and 
it lacks a user interface. 

• WHAMS - This code was developed by Professor Ted Belytschko of Northwestern 
University and represents a crashworthiness analysis program. The program features 
an h - adaptive procedure which automatically refines or coarsens the finite-element 
discretization in order to maintain solution accuracy. It uses subcycling and a contact- 
impact pinball penalty algorithm with an explicit time integration algorithm. 

• STAGS - This code was developed by Lockheed-Martin Palo Alto Research Labora- 
tory under NASA, Navy, and Air Force sponsorship. It represents a general-purpose, 
nonlinear shell, finite-element analysis code with static, transient dynamic, and eigen- 
value analysis options. The transient response is predicted using an implicit time 
integration procedure and only a static contact algorithm is available. If the impact 
force versus time curve is known, STAGS can be used to predict general trends and 
overall structural response. This code is available on NASA computers and through 
COSMIC. 

Perform Parametric Studies 

Parametric studies were performed in order to define the modeling and analysis needs 
for such simulations. These studies focused on simulating existing test data for low-speed 
impact problems using NASA experimental results and other results available in the lit- 
erature. Initially these studies were aimed at predicting the impact (or interface) force 
as a function of time assuming no damage to the laminate. The results presented here 
primarily represent a summary of the studies performed in Reference 5. All simulations 
were performed on a SGI Indigo II workstation with the R4400 processor. As a result, 
very long run times were experienced. 

The starting point for the impact simulations was a DYNA3D sample problem (Sample 
Problem 3) which is the impact of a cylindrical rod at the center of a square plate and the 
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supporting structure. The plate is isotropic with an elasto-plastic material model. The 
finite element model is shown in Figure 2 which assumes a doubly symmetric response 
from this quarter model. This model is referred to as the original model (or OM). A 
slightly modified form was considered wherein the physical supports were replaced by 
simple-support boundary conditions (or OM/SS). The baseline model (or BL) refers to 
the OM/SS model with a very high yield stress so that the response remains elastic. In 
this BL model, only one integration point through the thickness is used. Then a series 
of models denoted by LC i were examined where i increased from 2 to 7 indicating an 
increase in the number of integration points through the thickness. The structure is still 
isotropic so the results from BL and all LCi models should give the same response with 
the difference being the computational effort to obtain that response. 

The transient response of the transverse center deflection for each model is shown in 
Figure 3, The response obtained by changing support conditions is obvious and consistent 
with what was expected. The response for the other models indicate a vibration about 
the undeflected state after the impactor loses contact with the plate (i.e., rebound event). 
During the eaxly part of the transient, all models give essentially the same solution. How- 
ever, later on, differences are readily noted for models with few integration points through 
the thickness of the plate. As the number of integration points increases to 4, the transient 
response becomes more consistent. 

The computational effort for each of these simulations is given in Table 1. The column 
labeled Run Times is the wall clock time needed to perform the simulation in a non- 
dedicated mode (i.e., other jobs running at the same time possibly) and is not a clear 
measure of the increase in execution time as a result of increased model fidelity. For the 
LCi models, the same number of time steps is used. From this table, the increase in the 
number of integration points from 1 to 7 nearly doubled the required storage space on 
disk for the computational database and also nearly doubled the execution CPU time. As 
a consequence for laminated plates, the modeling of through-the-thickness effects using 
2-D plate/shell elements for a typical 8-ply laminate will result in a significant increase in 
computational effort unless ply clustering is used. 

The next study considered the same basic geometry for the plate and examined the 
effects of mesh refinement in the anticipated contact region. In addition, the impactor 
shape is changed to spherical rather than a flat end cylindrical rod. A schematic of this 
problem (Problem 3) is shown in Figure 4. Eight different finite element meshes were 
considered as indicated in Figure 5. Meshes 1 and 2 have uniform element spacings, while 
the other meshes have a graded mesh with more elements along the symmetry planes and 
in the contact region. Table 2 lists some general characteristics of each mesh. As the mesh 
is refined and element size decreases, the critical time step for the explicit time-integration 
method also decreases, and hence the total number of time steps needed to reach a specific 
simulation time (say 1000 micro-seconds) increases. This consequence is clearly indicated 
by the data given in Table 3. 

The transient response of the transverse deflection at the center of the plate and 
also at the point x = y = 2.5 inches from the plate center are shown in Figures 6 and 
7, respectively. The point directly under the impactor begins to move immediately (see 
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Figure 6), while the other point requires about 75 micro-seconds for motion to begin (see 
Figure 7). For the most part, all of the results from the different meshes exhibit the same 
overall trends, and results obtained using Mesh 7 are considered to be the more accurate 
results. 

The transient response of the surface strains on the lower surface of the plate at the 
point x = y — 2.5 inches from the plate center are shown in Figures 8 and 9. The strain in 
the x— direction is shown in Figure 8, and the strain in the y— direction is given in Figure 
9. Again the overall trends are similar, and the finite element mesh in the neighborhood 
of this point away from the center has approximately the same spatial discretization for 
all the graded mesh cases. 

The next simulation is for a laminated panel using the modeling strategies that have 
been developed from the earlier simulations of related geometries. A schematic diagram of 
this problem is shown in Figure 10, and this problem is referred to as Problem 4. The panel 
is a flat rectangular panel, and the impactor is a hemispherical dropped-weight assembly. 
The panel is a 48-ply quasi-isotropic graphite-epoxy laminate, and the impactor is steel. 
Material data are from Reference 8. Two different finite element models were used for 
this simulation as shown in Figure 11. The uniform mesh is considered to be the coarse 
mesh to reflect the modeling in the neighborhood of the impact, while the graded mesh is 
referred to as the refined mesh because of the high refinement in that same region. The 
smallest element size in the coarse mesh is ten times the smallest element size in the refined 
mesh. However, the total number of elements in the coarse mesh exceeds the number of 
elements in the refined mesh, while the number of elements in the contact region for the 
refined mesh is 25 times that in the coarse mesh. This has a direct consequence on the 
computational effort for the simulations (i.e., the simulation using the refined mesh will 
take approximately ten times as many time steps to reach the same amount of simulation 
time). 

The transient response for the predicted impact force versus time is shown in Figure 12. 
The results obtained from two DYNA3D simulations with different spatial discretizations 
are compared with the analytical results presented in Reference 8. The results obtained 
using the coarse mesh with DYNA3D are in reasonably good agreement with those of 
Reference 8 as indicated in Figure 12. The results obtained using the refined mesh ex- 
hibit a similar trend with increased amplitudes and a shorter impact/contact time. These 
differences may be attributed to the simulation fidelity of Reference 8, laminate model- 
ing, contact modeling, or a combination. Further studies will be needed to resolve these 
differences. 

The final simulation was performed using the LS-DYNA3D code. This simulation 
involved an oblique impact of a rod and a plate with both modeled only with solid elements. 
This problem formed the basis of a study to examine alternative ways to represent the 
contact surface and the associated events (e.g., rigid impactor). The results shown in 
Figure 13 are only qualitative and indicate the complexities of an impact/penetration 
simulation that will be needed for a composite structure. From Figure 13, the impactor is 
observed to deform and fragment as well as slide along the surface of the plate. In addition, 
the plate is punctured with through-the-thickness damage, and a hole is created by the 
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impactor, thereby requiring elements to collapse and/or erode. This type of capability is 
needed for simulation events such as the rotor burst problem. 

Identify Research Thrusts 

As a results of these studies, additional new research thrusts are defined. Simulation 
of the impact and penetration problem for composite structures involves many aspects of 
computational mechanics and composite mechanics that need to interact in a synergistic 
manner leading to an intelligent computational system for this class of problems. Issues as- 
sociate with damage growth, large deformations coupled with fragmentation, and residual 
strength predictions for composite wing-like structures include the need to: 

• Establish modeling criteria to develop crack simulation and growth models based on 
material response or fracture criteria. One such approach to developing a validated 
modeling strategy is included in Figure 14. 

• Identify failure modes and develop failure models for laminated composite structures 
including different failure mechanisms, fracture models, and material degradation 
models to be used in combination with point stress failure models. 

• Provide through-the-thickness damage modeling capability for laminated and sand- 
wich structures. 

• Provide an adaptive contact modeling algorithm that evolves as the impact and/or 
penetration progresses and contact surfaces change. 

• Provide robust and efficient procedure to establish the pre-stress state from which 
to initiate failure rather than simulate the penetration event. Adaptive dynamic 
relaxation is one procedure to evaluate. 

• Assess the structural response characteristics for assumed quasi-static stable crack 
growth simulations and those incorporating inertia effects from the impact and/or 
penetration event itself. Experimental tests on dynamic crack growth are needed in 
addition to numerical simulations. 

CONCLUDING REMARKS 

A preliminary investigation of simulating impact and penetration of composite struc- 
tures was performed under this grant. Existing finite element analysis codes were reviewed 
and their capabilities generally defined as related to the impact event. Literature review 
on the impact response of composite plates and composite curved panels was performed, 
and a detailed summary is reported in Reference 5. Limited detailed analysis simulations 
are available, and only limited experimental results are available. The analysis simulations 
of impact nearly always assume a contact force profile (or “footprint” of the impactor) as 
well as an impact force magnitude versus time distribution. In a design setting, neither 
of these attributes is generally available unless testing is performed. For flight hardware, 
such data are rarely available. 

Simulations performed to-date are limited to low-speed impact simulations which do 
not cause any material damage. These simulations indicate the modeling fidelity needed 
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to capture the transient dynamic response accurately, and these studies indicate the com- 
putational complexities to be anticipated for impact and penetration studies for laminated 
composites experiencing damage. 
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Table 1. Comparison of computational effort for Problem 2. 


Model 

Designation 

mS§ 

Number of 
Time Steps 

Disk Space 
(in Bytes) 


10 

Model Description^ 


| 0M (l> 

2.2 

12,191 


7,556 

26.62 

Original Model used as a BaseLine for comparison 


1.68 

12,163 

424,680 

6,155 

16.71 

Original Model with Simply Supported Boundary 
Conditions 

HflU 

3.45 

12,062 

424,680 

6,048 

44 84 

Original Model with No Supports and modified material 
properties 

LC2 ,!i 

4.38 

13,350 


7,747 

36.78 

Laminated Composite with 2 integration points 

LC3 {2> 

5.02 

13,350 


8,867 

44.69 

Laminated Composite with 3 integration points 

IQIH 

5.78 

13,350 

663,989 

10.165 

59.83 

Laminated Composite with 4 integration points 

LC5 <!> 

7.88 

13,350 

744,716 

11,735 

59.13 

Laminated Composite with 5 integration points 

LC6 (!1 

' 7.32 

13,350 

825,443 

12,844 

68.83 

Laminated Composite with 6 integration points 

Lcr 21 

3.82 

13,350 

906,170 

13,607 

38.32 

Laminated Composite with 7 integration points 


(1) Material properties as listed in Tabic 4 1 ? S' 

(2) Material properties as listed in Table 4.3 S 

(3) All simulation were run on the SGI in double precision 

(4) Data Dump Time interval was 1.05e+02 with 95 d3pk# files 

(5) Wall-clock time needed to run in a non-dedicated environment 
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Table 3. Comparison of computational effort for Problem 3. 


Mesh Name 

Smallest Ele. 

Run Time 

Number of 

CPU 

10 


mra mm 






Mesh 1 

0.2 

0.183 

2,301 

543 

8.05 

Mesh 2 

0.1 

3.5 


3,184 

25.03 

Mesh 3 

0.02 

10.7 

23,008 

16,016 

83.06 

Mesh 4 

0.01 

12.9 

46,016 

30,632 

107.18 

Mesh 5 

0.005 

20.5 

92,032 


212.04 

Mesh 6 ** 

0.00375 

65.8 

122,710 

114,022 

520.48 

Mesh 7 

0.0025 

62.5 


198,430 

756.68 

Mesh 8 *** 

0.002 

99.1 



1305.00 


* All simulations were run in a non-dedicated environment with the run time equal 
to the wall clock time. 

** This simulation was terminated at 5.0e-4 sec, therefore the computational results 
have been multiplied by 1.0e-3/5.0e-4 = 2.0 for scaling effects. 

*** This simulation was only run out to 9.0e-5 sec, therefore the computational 
results have been multiplied by 1.0e-3/9.0e05 = 1 1.11 for scaling effects. 


18 



































• Break tied nodes or coincident nodes concept 

similar to approach used for quasi-static analyses 
no mass loss 

assumed crack growth path 
break tie based on plastic strain 



• Eroded element concept 

requires small elements in the crack front 

needs transition modeling 

mass loss; function of element size 

“failed” elements removed or eroded from mesh 


Figure 1. Crack modeling concepts. 


Plat* Imp act ad by a rod (cm, ga, mtcrosao) 
tlm* - 0.00000E+00 








1 1 II! I 



x 1 

(b) Side View 


Figure 4. Panel dimensions for the quarter model of Problem 3. 
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Total Elements ■ 625 
Contact Region Element! - 4 

Contact Region Length divided by 
the Impactor Radius - ! 60 

Smallest Element Length * 0 2 in. 
Max Element Aspect Ratio - 1 



(a) Mesh t Shell Element Configuration and Characteristics 


Total Elements - 2.500 
Contact Region Elements * 9 

Contact Region Length divided by 
the Impactor Radius * 1 .20 

Smallest Element Length » 0. 1 in. 
Max Element Aspect Ratio * l 



(b) Mesh 2 Shell Element Configuration and Characteristics 



(f) Mesh 6 Shell Element Configuration and Characteristics 


Total Elements * 2,500 
Contact Region Elements * 25 

Contact Region Length divided by 
the Impactor Radius “ 0 40 

Smallest Elcm Length ■ 0.02 in. 
Max Ele. Aspect Ratio “ 5.S75 



(c) Mesh 3 Shell Element Configuration and Characteristics 


Total Elements - 4.225 
Contact Region Elements - 256 

Contact Region Length divided by 
the Impactor Radius “ 0.16 

Sm. Ekra. Length * 0.0025 in. 
Max Element Aspect Ratio * 60 



(g) Mesh 7 Shell Element Configuration and Characteristics 


Total Elements - 2.500 
Contact Region Elements - 25 

Contact Region Length divided by 
the Impactor Radius * 0.20 

Smallest Elem, Length “ 0,01 in. 
Max Element Aspect Ratio H 1 3 



Total Elements - 5.329 
Contact Region Elements - 225 

Contact Region Length divided by 
the Impactor Radius * 0. 12 

Smallest Elan. Length ■ 0.002 in. 
Max Element Aspect Ratio ■ 75 



(d) Mesh 4 Shell Element Configuration and Characteristics 


(h) Mesh l Shell Element Configuration and Cbartctermxa 


Figure 5. Mesh configurations used to simulate and study an impact event. 
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Transverse Deflection (inches) 




Deflection (in,) 



x-ii train (e^j (in./inj 



y-Strain (Eyy) (inyin.) 
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-2.00E-05 
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-8.00E-05 


— o— Mesh 1 

e Mesh 2 

— A — Mesh 3 

— g — Mesh 4 

Mesh if 

Mesh 6 

Mesh 7 

— I — Mesh 8 


SOT; 




Time (fi-scc) 


Figure 9. Transient response for the lower surface strain in the y-direction at the point 
located at x = y = 2.5 inches from panel center for Problem 3. 




v 0 = 88 in./sec 



(b) Side View 


Figure 10. Panel dimensions for the quarter model of panel impacted 
by a 2.6 lb m dropped weight assembly - Problem 4. 
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Figure 11. Mesh configurations used for Problem 4. 
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Total Impact Force (lb*) 



Figure 12. Impact force profile for Problem 4. 
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Figure 14. Crack growth modeling concept strategy. 
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